Purpose

QC metrics has been used to justify the pipelien quality. Here we propose an analysis to evaluate the pipeline updates/changes by using Picard QC metrics.

Input data

The key input data for this test are two Picardtools QC metrics from two pipelines

• bmetrics: QC metrics generated by base pipeline.
• umetrics: QC metrics generated by updated pipeline.
• metKeys: list of metrics key name to process
• output_name: output name string

First load functions we will use in this test

Here is the list of input parameters

source('~/Works/github/HCA/skylab/benchmarking/smartseq2/R/analysis_functions.R')
## inputs, python style
option_list <- list(
  make_option(
    "--bmetrics",
    type = "character",
    default = NULL,
    help = " First/Base metrics file name",
    metavar = "character"
  ),
  make_option(
    "--umetrics",
    type = "character",
    default = NULL,
    help = " Second/Updated metrics file name",
    metavar = "character"
  ),
  make_option(
    "--output_name",
    type = "character",
    default = "out",
    help = "output file name [default= %default]",
    metavar = "character"
  ),
  make_option(
    "--metKeys",
    type = "character",
    default = "out",
    help = "A list of metrics to analyze",
    metavar = "character"
  )
)
# arg rendered by rmkdown
args<-strsplit(commandArgs(trailingOnly = TRUE),split=' ')
opt_parser <- OptionParser(option_list = option_list)
opt <- parse_args(opt_parser,args=args[[1]])
# load params
output_name <- opt$output_name

Parse QC Metrics Files

The QC metrics files are the collection of Picardtools QC outputs, which are formated as N*M matrix The rows are QC metrics and columns are cells/samples.

met1 <- read.csv(opt$bmetrics, row.names = 1)
met2 <- read.csv(opt$umetrics, row.names = 1)

Here is the list of metrics we will analysis with

metKeys <- strsplit(opt$metKeys, split = ',')[[1]]
print(metKeys)

##  [1] "detected_ratio"               "MT_ratio"                    
##  [3] "PCT_PF_READS_ALIGNED"         "PCT_PF_READS_IMPROPER_PAIRS" 
##  [5] "PCT_READS_ALIGNED_IN_PAIRS"   "PCT_CODING_BASES"            
##  [7] "PCT_INTERGENIC_BASES"         "PCT_INTRONIC_BASES"          
##  [9] "PCT_UTR_BASES"                "PCT_USABLE_BASES"            
## [11] "PCT_MRNA_BASES"               "PCT_RIBOSOMAL_BASES"         
## [13] "PERCENT_DUPLICATION"          "MEDIAN_5PRIME_TO_3PRIME_BIAS"
## [15] "MEDIAN_3PRIME_BIAS"           "MEDIAN_5PRIME_BIAS"          
## [17] "MEDIAN_CV_COVERAGE"           "PF_MISMATCH_RATE"            
## [19] "MEDIAN_INSERT_SIZE"

First we match column names(cell IDs) between two pipelines metrics

colnames1 <- colnames(met1)
colnames2 <- colnames(met2)
mlist <- match(colnames1, colnames2)
nalist <- is.na(mlist)
if (sum(nalist) > 0) {
  print("input files have different column elements")
  exit()
}
met2 <- met2[, mlist]

Then select the subset of QC metrics to run analysis.

# select subset of metrics
met1.core <- met1[metKeys, ]
met2.core <- met2[metKeys, ]

Finally, we have two QC metrics to run analysis on.

Analysis

To evaluate the pipeline changes/updates, we will carry out statistical tests on each QC metric and these tests are described and visualized as following:

• Compare overall distribution in histogram.
• Run Wilcoxon test between two pipelines and visualize in violoin/boxplot. The pvalues are labeled. The significant p-value indicates the significant changes in this QC metric between two pipelines
• Linear regression test is carried out to examine the linearity of QC metric bwtween pipeline. The predicted model is presented in the xy-plot as well as the p-value.
• Kolmogorov–Smirnov(k-s) test is carried out to examine the similarity in distribution of QC metric between pipeline. K-S statistics(D-stats) is labeled as well as the p-value.

Let’s run the analysis mentioned above on the set of QC metrics.

• results=‘asis’, echo=FALSE

out <- c()
pouts <- list()
for (ii in 1:length(metKeys)) {
  x <- as.numeric(met1.core[metKeys[ii], ])
  y <- as.numeric(met2.core[metKeys[ii], ])
  z <- data.frame('Base' = x, 'Updated' = y)
  # linear regression model
  p1 <- plotlm(z)
  # plot histogram
  p2 <- plothist(z)
  # plot box
  p3 <- plotBox(z)
  # plot KS
  p4 <- plotKS(z)
  # test block
  ctext <-
    paste(
      "A: Density plot of",
      metKeys[ii],
      " from two pipeline.\n",
      "B: Violin of ",
      metKeys[ii],
      ", labeled with Wilcoxon test result.\n",
      "C: Linear regression test of ",
      metKeys[ii],
      "between two pipelines.\n",
      "D: K-S test, D statistics is labeled.\n",
      sep = " "
    )
  #cat(metKeys[ii],text,' \n')
  # arranage plots in one figure
  gp <-
    ggarrange(
      p2$p,
      p3$p,
      p1$p,
      p4$p,
      labels = c("A", "B", "C", "D")
    ) # arrange 4 plots
  gp <-
    gp + ggtitle(paste(metKeys[ii])) + theme(plot.title = element_text(
      hjust = 0.5,
      size = 15,
      face = 'bold'
    ))
  p<-annotate_figure(
    gp,
    bottom = text_grob(
      ctext,
      face = "italic",
      size = 13
    )
  )
  print(p)
  out <-
    rbind(out,
          c(
            metKeys[ii],
            p1$beta,
            p1$a,
            p1$r2,
            p1$fpval,
            p3$pval,
            p4$D,
            p4$pval
          )) # test, statistics outputs

  }

## Warning: package 'bindrcpp' was built under R version 3.4.4

Outputs

colnames(out) <-
  c('metrics',
    'beta',
    'a',
    'r2',
    'pvalue',
    'wilcoxon',
    'ks-D-stats',
    'ks-Pvalue')
write.table(
  out,
  file = paste(output_name, '_tests_stats.csv', sep = ''),
  quote = F,
  row.names = F,
  col.names = T,
  sep = ','
)